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ABSTRACT: Real-world systems are often complex, dynamic, and nonlinear. Understanding the dynamics 
of a system from its observed time series is key to the prediction and control of the system’s behavior. While 
most existing techniques tacitly assume some form of stationarity or continuity, abrupt changes, which are often 
due to external disturbances or sudden changes in the intrinsic dynamics, are common in time series. Structural 
breaks, which are time points at which the statistical patterns of a time series change, pose considerable chal¬ 
lenges to data analysis. Without identification of such break points, the same dynamic rule would be applied to 
the whole period of observation, whereas false identification of structural breaks may lead to overfitting. In this 
paper, we cast the problem of decomposing a time series into its trend and seasonal components as an optimiza¬ 
tion problem. This problem is ill-posed due to the arbitrariness in the number of parameters. To overcome this 
difficulty, we propose the addition of a penalty function (i.e., a regularization term) that accounts for the number 
of parameters. Our approach simultaneously identifies seasonality and trend without the need of iterations, and 
allows the reliable detection of structural breaks. The method is applied to recorded data on fish populations and 
sea surface temperature, where it detects structural breaks that would have been neglected otherwise. This sug¬ 
gests that our method can lead to a general approach for the monitoring, prediction, and prevention of structural 
changes in real systems. 


1 INTRODUCTION 


Systems in the real world are often complex not 
only with respect to the underlying interaction 
networks but also with respect to their dynam¬ 
ics. Examples include the temporal variations in 
economic growth (ICerra and Saxena 20081, fluc¬ 
tuations of metabolic rates dLabra et al. 20071) . 
oscillations in power-grid genera¬ 

tors (iFilatrella et al. 20081) . and dynamics of animal 
species (|Bj0rnstad and Grenfell 200T|). Understand¬ 


ing the dynamics from limited information (often in 
the form of a time series) is key for the prediction and 


control of system-level behavior. 

While time series analysis can benefit from mod¬ 
ern statistics, dynamical systems, and network theory, 
most existing techniques tacitly assume some form 
of stationarity or continuity (iBrockwell 20051) . How¬ 
ever, abrupt changes, which we refer to as structural 
breaks, are quite common in time series (Figure CO- 
These are often due to external disturbances or sudden 
changes in the intrinsic dynamics and pose consider¬ 
able challenges to data analysis and interpretation. 

Structural breaks are points in time at 
which the statistical patterns of a time series 
change (I Andrews 19931 IBai and Perron 20031) . With- 















out identification of such break points, the same 
dynamic rule would be applied to the whole period 
of observation, resulting in biases in the estimation 
of the system dynamics. On the other hand, false 
identification of structural breaks may split the time 
period into unnecessarily small subintervals, which 
affects statistical significance, introduces unnecessary 
parameters, and may lead to overfitting. 

In this paper, we cast the problem of decomposing 
a time series into its trend and seasonal components 
(traditionally achieved by an iterative scheme) as an 
optimization problem, whose objective function is the 
norm of the residuals. We show that this problem is 
ill-posed due to the arbitrariness in the number of pa¬ 
rameters. To overcome this difficulty, we propose the 
addition of a penalty function (i.e., a regularization 
term) that accounts for the number of parameters— 
a strategy often used in dealing with ill-conditioned 
linear systems (iNeumaier 19981) . This modest change 
leads to successful identification of seasonality and 
trend without the need of iterations. Furthermore, we 
show that our formulation allows the simultaneous 
decomposition of a time series into a trend compo¬ 
nent, a seasonal component, and a noise component, 
as well as structural changes in these components. 
Our approach is therefore a generalization of the clas¬ 
sical method of Bai and Perron (iBai and Perron 19981 
iBai and Perron 20031) . which only deals with time se¬ 
ries without a seasonal component, and is an al¬ 
ternative to the recently proposed methods for the 
identification of structural breaks in the trend of 


a seasonal time series (Haywood and Randal 2008 
iVerbesselt et al. 20l0 ). 

We validate our method using synthetic data and 
apply it to time series describing fish populations and 
sea surface temperature. We found structural breaks 
that would have been neglected using previous meth¬ 
ods. This indicates that our method is promising in the 
development of improved approaches for the monitor¬ 
ing, prediction, and prevention of structural changes 
in real systems. 



time, t 

Figure 1: Example of a noisy time series that exhibits structural 
changes in both the trend component (dashed lines) and seasonal 
component (arrow). 


2 GENERAL STRUCTURAL BREAK MODEL 

The temporal dynamics of a system is often repre¬ 
sented by a time series {Y t }f =1 , where Y t e M is 
the state of the observed variable at time t. The 
classical seasonal-trend decomposition of a time se¬ 
ries Y t can be expressed as (iCleveland et al. 19901 
iBrockwell 20051) 

Y = 7t + S t + S t , (1) 

where % is the trend component (usually modeled as 
a function of t), S t is the seasonal component ( S t+d = 
S t where d > 0 is the period of the component), and 
St is random noise (which is assumed to have zero 
mean). 

Although Eq. (OQ) is generally valid, finding the de¬ 
composition itself turns out to be a nontrivial task. In 
particular, in the presence of structural breaks the de¬ 
composition will depend on the location of the break 
points, which are themselves dependent on the de¬ 
composition. This intrinsic coupling renders classi¬ 
cal decomposition techniques (ICleveland et al. 19901 
iBrockwell 20051) inappropriate in general. 

To incorporate the presence of possible structural 
breaks into the trend component, we divide the time 
interval [0, T\ into m subintervals according to the 
partition 

0 = t*<f*< ...<4 = T, (2) 

and assume that the trend (as a function of time) re¬ 
mains the same within each subinterval t* + 1 < t < 
t *_|_i (i — 0,1,... ,m — 1). Similarly, for the seasonal 
component, we partition the time period into n subin¬ 
tervals according to 

0 = t+ < t+ < • • • < t+ = T, (3) 

and assume that the seasonal component is un¬ 
changed within each subinterval tf + 1 < t < t+ +1 
(i = 0,l,...,ra-l). 

Here we model the trend in the Ath trend subinter¬ 
val as a linear function of time, 

71 = (i t t + bi, t* + 1 <■ t < (4) 

On the other hand, the seasonal component in the A 
th seasonal subinterval is represented by a set of di 
numbers {s^, s%\ • • •, s^}, as 

= tt + l<t<t+ +1 , (5) 

with k = t — di[t/dj\, where di is the period of sea¬ 
sonality in the /-1h subinterval and |ycj is defined as 
the largest integer that is smaller than or equal to x. 
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3 SEASONAL-TREND DECOMPOSITION 


We first consider the problem of seasonal-trend de¬ 
composition of a time series {Y t }J =1 in the absence 
of structural breaks. In particular, the goal is to de¬ 
compose {Y,}J =i as in Eq. ([Q) and at the same time 
estimate the parameters a and b in Eq. © for the 
trend component and the parameters {si, s 2 ,..., s d } 
in Eq. © for the seasonal component. 

3.1 Traditional Approach 

The traditional approach for seasonal-trend decompo¬ 
sition often involves three steps (IBrockwell 2005b . 

The first step is to estimate the trend % by applying 
a moving average filter of size q that attempts to elim¬ 
inate the seasonal component and meanwhile reduces 
noise. The resulting series is computed as 

1 t+[q/2\ 

T> = - E W ’ Y ’• < 6 > 

i=*-L<j/ 2 J 

where the coefficients are w = [1,1,..., 1,1] if q is 
odd and w = [0.5,1,..., 1,0.5] if q is even. If the pe¬ 
riod d is given, a natural choice would be q = d. 

The second step is to estimate the seasonal compo¬ 
nent. For each k — 1,2,..., q, we compute s fc as 


s k (Yk+jq) ^ 


(7) 


where (•) indicates average, and j takes all integer val¬ 
ues for which 1 < k + jq < T. Once the values of Sk 
are calculated, we obtain the seasonal component S t 
according to Eq. © for the entire time period. 

The third step is to re-estimate the trend compo¬ 
nent, often by fitting a least squares polynomial to the 
series {Y t — S t }. The resulting polynomial is taken as 
the trend component Tt, and its least squares polyno¬ 
mial fit gives the model parameters used in Eq. ©. 

In Figure [2] we show that the above approach is ef¬ 
fective when the size of the moving average filter is 
chosen to be q = d, and may generate misleading re¬ 
sults when a different value of q is used. Since d is 
unknown in general, this approach often requires trial 
and error before an appropriate decomposition can 
be achieved. We note that there are other methods to 
perform seasonal-trend decompositions, including the 
STL procedure proposed by [Cleveland et al. (1990] ). 
Most of the more sophisticated methods require iter¬ 
ations of multiple steps, and yet do not guarantee a 
reliable decomposition into seasonal and trend com¬ 
ponents. 



Figure 2: Classical seasonal-trend decomposition, (a) Syn¬ 
thetic time series generated by the equation Y t = at + b + 
sin(27 xt/d) + e t , where a = 0.03, b = —0.5, d = 10, and each 
£t is generated independently from a Gaussian distribution with 
mean p = 0 and standard deviation a = 0.1. (b) Actual seasonal 
component^ = sin(27rf/d) and its estimated counterpart by the 
traditional approach described in Section l3TTl with moving aver¬ 
age biter size q = d in Eq. ©. (c) Same as in panel (b), but now 
for q = 1.5 d. 


3.2 Regularized Optimization Approach 

We propose an optimization-based approach for the 
seasonal-trend decomposition. The idea is that, for a 
given trend model (e.g., the linear model in Eq. ©) 
and estimated seasonal period p, we have the follow¬ 
ing matrix equation for Y = [>'i, Y 2 ,.... Y T ] J : 

Y = Q5 + S 


" 1 1 1 ... O' 





a 



b 

p 1 0 ... 1 

p + 1 1 1 ... 0 


Sl 


S2 

2 p 1 0 ... 1 


Sp 
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where Q is a T x (p + 2) constant matrix, 5 is a (p + 
2) x 1 vector of parameters, and £ is a T x 1 residual 
vector. The problem of decomposing the time series 
into its trend and seasonal components then becomes 
the problem of estimating the parameter vector S in 
Eq. ([8]). A common criterion is to minimize the sum 
of the squares of the residuals, ||£j| 2 . For fixed p and 
Euclidean norm, || • || 2 , this criterion leads to the least 
squares estimate: 

5 = Q + Y, (9) 

where Q + is the pseudo-inverse of the matrix 
Q (IHorn and Johnson 19851) . 

Note that, if we are allowed to freely choose p, the 
choice p = T and § = [0, c, Y J — c] T would yield £ = 
0 for any number c, which corresponds to the minimal 
possible sum of squared residuals. In other words, the 
problem of minimizing ||£|| has an infinite number of 
solutions and is therefore ill-posed as is. 

To overcome this difficulty, we propose the addi¬ 
tion of a penalty function to the objective function. 
For the given Y, we solve the regularized optimiza¬ 
tion problem 

min J(5) = ||F — QS\\ 2 + A||5|| 0 , (10) 

5eKP+ 2 , 0 <p<T 

where Q is the matrix defined in Eq. ([HI), A > 0 is 
a (predefined) regularization parameter, and A||<5|| 0 is 
the penalty term where ||£|| 0 accounts for the number 
of parameters in the model. 

For a given regularization parameter A, the opti¬ 
mization problem (fTOl) can be solved in two steps: 

Step 1: For each p, find the solution 5^ 

that minimizes J(5) for A = 0. ,,,, 

v ' (11) 

Step 2: Choose 5 = argmin 0<p<T J(<5 (p )) 

for the given A. 

In practice, the regularization parameter A is often 
chosen to be a positive number that is small relative 
to the sum of squared residuals term in Eq. (fTOl) . 

We apply this regularized optimization approach to 
the synthetic time series used in Fig. [2] Figure [3j a) 
shows that the regularized optimization successfully 
detects the true period of the seasonal component 
when A = 0.1. In fact, the minimal value of J is 
achieved at p = d = 10 for any A G (A min , A max ), 
where, in this example, A min ~ 0.01 and A max ~ 0.5. 
This suggests robustness of the regularized optimiza¬ 
tion approach with respect to the regularization pa¬ 
rameter A. Figure [3lb) shows the excellent agreement 
between the estimated seasonal component and the 
actual data (the fit to the trend, which is not plotted, is 
also excellent for this example). 




Figure 3: Seasonal-trend decomposition by regularized opti¬ 
mization. The time series is the same shown in Fig. [ 2|a). (a) 
Dependence of the objective function J(S W) on p, where J is 
defined in Eq. ( ITOl ) for A = 0.1 and 6 tj,) denotes the optimal 
solution for the given p. Function J attains its minimum when 
p = d = 10, which is the true period of the seasonal component 
in this example, (b) Actual seasonal component St = sin(27rf/d) 
and its estimated counterpart by the regularized optimization ap¬ 
proach summarized in Eqs. (10-11). 

4 STRUCTURAF BREAKS IN TREND 

We now turn to a slightly different problem, where 
the goal is to identify structural breaks in the trend 
component in the absence of seasonality. In particular, 
we assume that S t = 0 in Eq. 0, but m > 1 in Eq. ©. 
In this case, the model for Y becomes 

Y = Q6 + £ (12) 


T t* 0 +1 





do 



b 0 

1 t\ 


CL\ 



h 

1 C-1 + 1 





Q"m—1 

1 c . 


1 _ 


where Q is T x 2 m and 5 is 2 m x 1. For a given time 
series {Y}J =1 , the ultimate goal within the modeling 
framework of Eq. (fl2l) is to estimate the number of 
segments m, the break points and the corre¬ 

sponding parameters ( ak,bk ) within each segment. 
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4.1 Dynamic Programming Approach 


4.3 Example: Fish Populations in Green Bay 


Model (fl2l) has the same form as the pure struc¬ 
tural change model (IBai and Perron 19981) . For a 
fixed number of segments m, an efficient way to 
find the break points and the corresponding param¬ 
eters that minimize the sum of residual squares is 
to use dynamic programming (iBai and Perron 20031) . 
Let SSR (i,j) be the sum of squared residuals ob¬ 
tained by applying least-squares to a segment that 
starts at t — i and ends at t = j, i.e., 

3 

SSR (i,j) = min (Y t — a — bt) 2 . (13) 

a.b ' ^ 
t—i 

Furthermore, let SSR({r; k}) be the minimum sum 
of squared residuals for the first r values of Y using 
k breaks. The desired solution SSR({T; m}) satisfies 
the following recursive equation 

SSR({T;m}) = 


min [SSR({j ;m — 1}) + SSR(j + 1, T)}. (14) 

1 <j<T 

Note that in the calculation of SSRone 
can use the recursive relation that relates 
SSR (i,j) to SSR (i,j — 1) for computational ef¬ 
ficiency (IBrown et al. 19751) . In a time series of 
length T, the dynamic programming algorithm 
involves order T 2 operations. Although not to be 
discussed in this paper, we point out that there has 
been recent work addressing in detail computational 
aspects of the dynamic programming algorithm and 
how its efficiency can be improved (Rigaill 2010). 


4.2 Choosing the Number of Breaks via 
Regularized Optimization 

The dynamic programming approach requires the 
number of breaks m, which in general cannot be 
determined a priori. Instead of relying on statisti¬ 
cal tests (IZeileis et al. 20031 IZeileis 2005!) . which of¬ 
ten assume a specific form for the distribution of 
the residual series, here again we treat the problem 
(of determining the number and location of structural 
breaks) as a regularized optimization problem: 

minimize ( TSSR({T;m})l + 2m \], (15) 

{t* k y£ =1 ,l<m<T/2\ L J J 

where the SSR term accounts for the sum of squared 
residuals, A > 0 is a regularization parameter, and 2m 
accounts for the total number of parameters when the 
time series is partitioned into m segments. 


We apply the regularized optimization approach to the 
population abundances of 43 fish species in the Green 
Bay, Wisconsin. The time series for each species 
contains 30 data points, corresponding to the annual 
abundance of that species from 1980 through 2009. 

Figure @Ja) shows the population abundance of 
the fish sheepshead and its best linear trend. Using 
A = 0.15 in Eq. (fl5l) . we find that the optimal so¬ 
lution requires m — 2 segments, with the structural 
break occurring at t\ = 12 (year 1991), as shown in 
Fig. [4(b). Furthermore, we apply Eq. (fl5l) to all 43 
species with A = 0.15 and obtain a distribution of the 
break times for all species as shown in Fig. [4(c). It 
is interesting to observe that the distribution curve 
has two peaks: one in the early 1990’s which is the 
time the invasive species round goby was first dis¬ 
covered in the system, and another peak in the early 
2000’s (iLederer et al. 20061 iLedereret al. 20081) . 



Figure 4: Structural breaks in the Green Bay fish populations, 
(a) Annual population abundance of sheepshead, a fish species 
from Green Bay. The solid line shows the least squares fit of the 
time series, (b) Same time series as in (a), partitioned into two 
segments using the dynamic programming algorithm described 
in Section 14.11 The vertical dashed line marks the break point 
and the two solid lines represent the least squares fit for the two 
time series segments, respectively, (c) Probability distribution of 
break points for all 43 species identified using our method. Note 
the peaks around 1993 and 2001 of this distribution. 
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5 STRUCTURAL BREAKS IN THE TREND OF 
A SEASONAL TIME SERIES 


as well as the corresponding parameters 5r and 5 S de¬ 
fined in Eqs. (fl6l)^([l9l): 


We now present our approach to the problem of de¬ 
tecting structural breaks in a time series that contains 
both trend and seasonal components. Such a time se¬ 
ries can be decomposed using the general form in 
Eq. ([]]) with the possible presence of structural breaks 
defined by Eqs. ([2}|5]). 


5.1 Regularized Optimization Formulation 

Using the notation introduced in Eqs. (Q]j5]) for the 
general seasonal time series with structural breaks, 
the time series model for Y can be expressed as 


Y = Qt$t + Qs$s + £■ 
Here, the term 

T t* + 1 


Qt$t — 




1 C-1 + 1 

1 t* 


(16) 



a i 
bi 


(17) 


1 

.bm— 1 . 


denotes a linear model for the trend component with 
m — 1 structural breaks (m segments). The term 


Qs$s 




Q 


(n) 



r s (1) i 


s ( 2 ) 

J Tx Ed, 

_s (n) _ 


(18) 


E d .. xi 


models the seasonal component with n — 1 breaks (n 
segments), where 




( 20 ) 


where Q and 5 are defined as 

Q = [Qt,Qs\, 3 = [dr 1 , £s T ] T - (21) 

Although the global optimum of problem (l20l) can 
be found by enumeration of all possible break points 
by solving a least squares problem for each set of 
break points, such a brute-force approach is not prac¬ 
tical due to its high computational cost. Here we focus 
on an alternative method, which is based on iterative 
optimization of the trend and seasonal parameters, 
respectively. In particular, we consider the scenario 
where the seasonal period d does not change, but al¬ 
lows for the presence of seasonality, trend, and struc¬ 
tural breaks in the trend component. Under such con¬ 
ditions, the following procedure is adopted for the de¬ 
tection of structural breaks and estimation of parame¬ 
ters (including those for the seasonal component): 


" Step 0: Initial assignment of T =0. 

Step 1: Seasonal-trend decomposition of Y — T 
via the methods from Sec. !3.2l to obtain S. 

Step 2: Estimation of the structural breaks and 
parameters of Y — S using the methods 
described in Secs. 14.11 - 14.21 
Update the estimated trend as T — Qt ( W- 

( 22 ) 

Steps 1 and 2 are repeated until the estimation of the 
breaks points and trend converges. 


Q 




T 

... O' 

0 

... 1 

1 

... 0 

0 

... 1 




' (*)' 
b 2 


- S d 


(0 


E d , xi 


df-tf- i)xdi 


(19) 


is the seasonal model for the ith segment (i = 
1)2,... ,n). 

For the given Y, and regularization parameter A > 
0, our approach is based on solving the following reg¬ 
ularized optimization problem to find the break points 


5.2 Example: Arctic Sea Surface Temperature 

We apply our method to the time series of the sea sur¬ 
face temperature (SST) from the Arctic. The time se¬ 
ries data is obtained from the U.S. Geological Sur¬ 
vey website (USGS 2013 ). The original data contains 
moderate-resolution monthly SST from 20 regional 
seas in the Arctic over 28 years (1982-2009). 

We focus on the overall monthly SST of the Arctic 
obtained by averaging over the SST from all 20 re¬ 
gional seas. The data for several regional seas are not 
available for the entire time period. To address this 
issue, we fill in the missing data by the temporal av¬ 
erage SST from available data for that sea. After this 
pre-processing step, we obtain the time series {Y t }J =1 , 
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(2002-2009) is decreasing rather than increasing. 
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Figure 5: Structural break in the Arctic sea surface temperature 
time series, (a) Monthly SST from January 1983 through De¬ 
cember 2009. (b) Seasonal-adjusted time series, Y — S, as well 
as its linear dependence of time (solid straight line), which re¬ 
sembles the global trend of the time series, (c) Same as (b), 
except that a break point around 2003 is identified using our 
approach (vertical dashed line). The break point separates the 
seasonal-adjusted time series into two temporal regimes, one 
slowly increasing in time and the other slowly decreasing in time 
(solid straight lines on the two sides of the breakpoint). 


as shown in Fig. [51a), where Y, denotes the average 
SST of the Arctic during the t-th month (t — 1 corre¬ 
sponds to January 1983 and T = 336). 

Using a regularization parameter A = 0.1, we found 
that the optimal seasonal period for the time series 
is d = 12, with the presence of 1 break point in the 
trend component, around the month of April in year 
2003. Figures |51b-c) show the seasonal-adjusted time 
series as well as their trend component when there 
is no structural break and when there is one struc¬ 
tural break (found to yield optimal solution via our ap¬ 
proach). Note that the current trend of the Arctic SST 
would have been asserted to be increasing if structural 
break in the trend component were ignored. However, 
our result suggests that, despite the overall increase 
for 20 years (1982-2001), the recent trend of the SST 


6 CONCLUSIONS 


In this paper, we proposed a regularized optimiza¬ 
tion framework for time series analysis, including re¬ 
newed approaches for solving the classical problem 
of seasonal-trend decomposition and the detection of 
structural breaks in a time series with or without the 
presence of seasonal components. Our approach was 
tested against synthetic data and applied to empiri¬ 
cal time series, including fish populations from Green 
Bay and sea surface temperature in the Arctic. We are 
able to detect structural breaks that are of practical 
importance in the prediction of future states of these 
complex systems. 

The key in our approach is the formulation that 
naturally treats parameters for the trend and sea¬ 
sonal components in a similar manner. The regular¬ 
ization term is analogous to the Akaike information 
criterion (IAkaike 19741) and the Schwarz-Bayesian 
criterion (ISchwarz 19781) . both of which are popu¬ 
lar methods for model selection. In our formula¬ 
tion, these methods can be useful for the selection 
of regularization parameters. Other approaches, such 
as the detrended fluctuation analysis (IHu et al. 20011 
IChen et al. 20021) . are also likely to be useful when in¬ 
tegrated into our framework. 

We used Euclidean norm to measure the model 
quality, due to its analytical convenience (solution 
to least square problems can be found easily us¬ 
ing singular value decomposition). In view of the 
recent developments in parameter estimation un¬ 
der general p norms, including p < 1 for sparse 
data ( Candes et al. 2006 ). and p = oo when the under¬ 
lying dynamics is chaotic (ISun et al. 2011) . it will be 
interesting to explore such scenarios and how they can 
improve the estimation of break points. Faster compu¬ 
tational procedures, however, will be necessary since 
most such problems involve much higher computa¬ 
tional cost than ordinary least square problems. 

Finally, the problem of predicting future occur¬ 
rences of structural breaks—or tipping points, in 
the language of climatology and ecology—remains a 
challenging problem (Scheffer et al. 20091) . It is our 
hope that the development of new and improved 
methods for the prediction of future trends, seasonal 
patterns, and breaks will benefit from from methods 
(such as the one proposed here) to estimate those 
quantities reliably in recorded time series. 
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